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We have characterized the structure of 176 different single-layer graphene grain boundaries using 
>1000 experimental HRTEM images using a semi-automated structure processing routine. We 
introduce a new algorithm for generating grain boundary structures for a class of hexagonal 2D 
materials and use this algorithm and molecular dynamics to simulate the structure of >79 000 
graphene grain boundaries covering 4122 unique orientations distributed over the entire parameter 
space. The dislocation content and structural properties are extracted from all experimental and 
simulated boundaries, and various trends are explored. We find excellent agreement between the 
simulated and experimentally observed grain boundaries. Our analysis demonstrates the power of 
a statistically significant number of measurements as opposed to a small number of observations in 
atomic science. All experimental and simulated boundary structures are available online. 


Introduction 

Single-layer graphene is a promising material for many 
technological applications, due to its excellent mechani¬ 
cal [1-3] and electronic properties [4, 5]. Most graphene 
deposition methods produce polycrystalline sheets, con¬ 
taining grain boundaries (GBs) [4, 6]. This polycrys¬ 
tal structure introduces local property deviations at the 
boundaries that could limit or enable various potential 
applications. There is also strong scientific interest in 
graphene GBs due to their one-dimensional nature. Some 
examples include a bimodal phonon scattering behaviour 
across graphene GBs [7] , anomalous strength characteris¬ 
tics [3, 8], strong chemical sensitivity of boundary charge 
transform [9], a transformation of the GBs from trans¬ 
parency of charge carriers to near-perfect reflection [10], 
amongst others. 

A large number of theoretical studies on graphene GB 
structures have been performed by many researchers [10- 
30]. However the number of corresponding experimental 
measurements of free-standing graphene GB structure at 
atomic resolution is much smaller [3, 31-35]. These ex¬ 
perimental studies have typically considered either a sin¬ 
gle boundary structure, or a small number of GB struc¬ 
tures. Thus, the gap between the number of possible 
or proposed graphene GB structures and those experi¬ 
mentally observed is very large. This makes testing the 
various proposed models for graphene GB structure and 
structure formation very challenging [11, 36]. 

In this study, we have characterized the structure «176 


graphene GB structures from 1067 phase-contrast high 
resolution transmission electron microscopy (HRTEM) 
observations of free-standing single-layer graphene sam¬ 
ples. We have characterized the atomic structure using 
a semi-automated processing routine, measuring the lo¬ 
cal topology and various other physical parameters. We 
have also used a new algorithm to generate the structure 
of ~79,000 graphene GBs covering the entire orientation 
parameter space, which were then relaxed using molecu¬ 
lar dynamics (MD). We have performed a detailed struc¬ 
tural characterization of all experimental and simulated 
boundaries, extracting structure parameters and disloca¬ 
tion content of all boundaries. The proposed algorithm 
for generating graphene GB structures is found to be in 
excellent agreement with the observed structures. 

Methods: Experimental 

Graphene Sample Fabrication and HRTEM Imaging 

Graphene samples were grown on polycrystalline cop¬ 
per substrates at 1035°G by chemical vapor deposition. 
The substrate was held at 150 mTorr hydrogen for 1.5 
hours in a quartz tube furnace, and then 400 mTorr 
methane was run (5 seem) over the sample for 15 min¬ 
utes to grow single-layer graphene. Further details of this 
method are described in Refs. [3, 4, 6]. 

All experimental high-resolution transmission elec¬ 
tron microscope (HRTEM) images used in this study 
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were recorded on the TEAM 0.5, a monochromated, 
aberration-corrected FEI Titan-class microscope, oper¬ 
ated at 80 kV with a high brightness Schottky field emis¬ 
sion gun. Rather than optimizing the imaging condi¬ 
tions, we instead focused on recording images as quickly 
as possible so as to minimize electron beam damage of 
the GBs. Often, multiple images of the same boundary 
were collected sequentially which allowed for both opti¬ 
mization of the imaging conditions and observations of 
beam-induced modifications of the structure. 

Semi-Automated Experimental GB Analysis 

The graphene HRTEM images have a low signal-to- 
noise ratio due to the low scattering efficiency of single 
carbon atoms. In order to measure the boundary struc¬ 
ture for hundreds of images, we have developed a pro¬ 
cessing routine with as few user inputs as possible. This 
routine is shown schematically in Fig. 1. First, a lin¬ 
ear background is fitted and removed from the image to 
minimize the intensity falloff caused by the monochro¬ 
matic aperture. Next, a nonlinear filter is applied to 
remove noise (rank filter of local intensity values within 
a O.bA radius, 25% darkest value selected [37]), shown 
in Fig. IB. Peak detection is used to find local minima, 
and the user inputs the boundary extent, the three loca¬ 
tions labeled in Fig. 1C. A subset of the detected peaks 
given by this user-selected region-of-interest is used to 
automatically characterize the boundary structure. 

The first step of the boundary characterization is a 
Voronoi tessellation of the detected local minima, shown 
in Fig. ID. The Voronoi cell vertices represent atomic po¬ 
sitions. The number of carbon atoms in each ring is given 
by the number of sides of each cell. Next, the boundary 
cells are removed and the tessellation is recomputed using 
a weighted Voronoi algorithm [38] , with the weights set to 
keep the mean bond length constant for all cells, shown 
in Fig. IE. The final atomic coordinates are plotted in 
Fig. IF, with a final optimization performed to enforce 
a minimum distance constraint of 1.2 A on all atomic 
coordinates, to ensure a physically realistic result. The 
boundary can be traced by connecting all pentagon and 
heptagon rings, and a best-fit lattice is calculated for the 
two grains on each side. The dislocations are located 
by searching for a minimal description of all pentagon- 
heptagon pairs. 

Additionally, the strain of the experimental boundaries 
was estimated using a geometric method similar to that 
given in Ref. [39]. Each atom is placed at the center 
of a triangle defined by its 3 nearest-neighbors, calcu¬ 
lated from a Delaunay triangulation of the set of atoms. 
These triangles are referenced to the appropriate trian¬ 
gle (2 atomic sites per unit cell) formed by the lattice 
vectors of the best-fit lattices of the two grains. The 
linear transformation matrix for each triangle is used to 



FIG. 1. Example of GB structure parsing for an experimental 
image: (A) original micrograph, (B) local intensity-ordered 
filtered image, (C) ring center detection and user-selected 
boundaries, (D) Voronoi tessellation, (E) edge cleanup and 
atomic coordinate generation, and (F) final coordinates with 
dislocation identification. 


calculate local strains (infinitesimal strain is assumed for 
convenience). Rather than parse the strain into atomic 
coordinates as in Ref. [39] , we instead calculate the root- 
mean-square values of the local strains and the local ro¬ 
tation, since we are interested in the “average” distortion 
of each of the boundaries. Three examples of these strain 
measurements are given in Fig. 2. 

Methods: Numerical 
Space of Graphene GBs 

Bulk three-dimensional materials require 5 angles to 
characterize the macroscopic degrees of freedom of a gen¬ 
eral GB, while two-dimensional materials require only 
2 angles. Thus the parameter space for 2D GBs is far 
smaller than that of 3D grain boundaries. As shown in 
Fig. 3 A these two angles are the misorientation angle 
+ ^2 5 defined as the angle between the unit cell 
vectors of each grain, and the boundary line direction 
= 1^1 — ^ 2 1, defined as the angle between the bound¬ 
ary vector and the symmetric tilt boundary vector for a 
given ^M- Due to the symmetries of the graphene lattice 
we get Om ^ (0°,60°) and 6 >l G (0°,^m)- A third pa- 
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FIG. 2. Examples of atomic strain measurements on three 
experimentally measured boundaries. The values of Cxx, Cyy, 
and Cxy range from -25% (blue) to +25% (red), while local 
rotation 0 has a range of ±20° from blue to red. 



FIG. 3. (Golor online) (A) A typical grain boundary struc¬ 
ture showing the graphene lattice vectors vi and V 2 , and the 
definition for the grain angles (^i, ^ 2 ). The misorientation and 
the line angles are defined as +^ 2 , and = \ 0 i— 02 \. 

(B) A representation of all the (^m,^l) pairs simulated in 
this study. The red dots indicate the perfectly commensu¬ 
rate GBs, while the black dots indicate the approximately 
commensurate ones. A total of 4122 (^m,^l) pairs were sim¬ 
ulated, resulting in over 79,000 GB structures when including 
the displacement A. 


rameter, namely the relatively sliding of the two grains 
along the GB boundary is also needed for a complete de¬ 
scription of the boundary. In our simulations we choose 
the relative sliding that gives the lowest GB energy, thus 
effectively eliminating this degree of freedom. 

In order to minimize the boundary effects, we simu¬ 
late GB structures that are periodic along the GB direc¬ 
tion. This requirement places strong restrictions on the 
GB configurations that we can simulate. Consider sim¬ 
ulation of the GB corresponding to a point or 

equivalently (^ 1 , O 2 ), in the parameter space. The lattice 
vectors for graphene are vi = a\/3ei, V 2 = a\/3/2ei + 


3a/2e2, where a is the carbon-carbon bond length, and 
ei ^2 are unit vectors parallel and perpendicular to the 
zigzag axis of the graphene sheet, respectively. Thus for 
given (^ 1 , O 2 ) the corresponding grains have periodic re¬ 
peat distances of li = a^/3{n‘f^ + n ?2 ± (^ = 2) 

along the GB direction, where 71^2 are integers such 
that tan6>i = (2nii +n^2)/\/3^i2 [40]. For an arbitrary Oi 
there may exist no suitable integers 77 ^ 2 , or even if 
such integers exist, li can be prohibitively large for MD 
simulation. Further, in order to simulate a GB, li and I 2 
should either be commensurate, i.e., I 1 /I 2 = p/q where 
p, q are positive integers, in which case the net GB length 
is given hy Iqb = Qh = ph^ or they should be approx¬ 
imately commensurate, i.e., h/h ~ p/q, in which case 
the simulated GB length is Iqb = ‘^hQhp/{hQ-\-hp)^ and 
each grain has a strain of magnitude Ihq — hpl/ihQ-i-hp)- 
In case of approximately commensurate boundaries, we 
require that the rational approximation p/q is chosen 
such that the resulting strain magnitude is less than 10“^, 
so that the resulting elastic distortion is minimal. Due to 
numerical considerations, we simulated GBs with a max¬ 
imum length of 2000 A. If for a given (6>i, O 2 ) the GB 
length is greater than 2000 A, we try to find a nearby GB 
such that the resulting grain angles are within 0.01° of 
the desired values. Finally, in order to choose the relative 
sliding between the two grains that leads to minimum GB 
energy, we search in steps of 1 A over the entire range 
A [41], given by 


( h 


A = 


min (|^2 - Ih/hlhl, 

Ih — R2Ai1^i|), 


if commensurate 
if approximately 

commensurate 


( 1 ) 


assuming li < I 2 . We simulate all perfectly commen¬ 
surate GBs with length less than 2000 A , and grid 
the (6>i, O 2 ) space in steps of 0.5°, resulting in a ‘di¬ 
agonal gridding’ of the (6 >m, ^l) space in steps of 1.0°. 
However, for certain configurations near (6 >m, ^l) = 
(0°, 0°), (60°, 0°), (60°, 60°) no approximately com¬ 
mensurate boundaries with length less than 2000 A could 
be found, and thus no boundaries were simulated at these 
grid points. Fig. 3B shows (6 >m, ^l) for all GBs config¬ 
urations that we simulated (4122 total). Each point in 
that figure represents several simulations due to the sam¬ 
pling of the relative sliding A. In all we have simulated 
over 79,000 GB structures for this study. 


Numerical GB Structure Generation Algorithm 

Experimental observations of well annealed graphene 
GBs in the present study, as well as by several previ¬ 
ous authors, [3, 31-35, 42, 43] suggests that the orien¬ 
tation change between adjacent grains meeting at a GB 
is mediated largely by pairs of rings of 5 and 7 carbon 
atoms. These pent agon-heptagon pairs, also called the 5- 
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7 pairs, are the dislocation cores with the shortest Burg¬ 
ers vectors in graphene, and have a low formation en¬ 
ergy [11]. Thus, it is reasonable to expect that graphene 
GBs simulated on a computer have atomic structures 
where the orientation change between the grains is me¬ 
diated mostly by the experimentally observed pentagon- 
heptagon pairs. However, it is difficult to meet this re¬ 
quirement in practice. While a few simple GBs com¬ 
posed solely of pentagon-heptagon pairs have been simu¬ 
lated successfully [21-25], deviations from this motif are 
evident in the GB and polycrystals used in several re¬ 
cent studies [26-30]. The reason for this limitation is 
that so far no computationally efficient method has been 
proposed to generate well-annealed graphene GBs on a 
computer. Methods based on grand canonical Monte 
Garlo simulations, while theoretically sound and simple 
to implement, take inordinately large amount of com¬ 
puter time in practice. The common method of simply 
“annealing” a grain boundary by running molecular dy¬ 
namics at an elevated temperature is also not effective 
since the typical thermal barriers that need to be over¬ 
come for suitable reconstruction are high, and thus the 
desired annealing does not occur in the limited simu¬ 
lation time. For simple GBs these limitations can be 
overcome by the inclination-disclination based geometric 
method [11]. This geometrical approach is well suited 
for the study of simple GBs, but becomes unwieldy for 
tailoring more complex GBs or polycrystals with several 
different GBs and junctions. 


FIG. 4. (Color online) (A) A 2D-periodic Voronoi tessellation 
with n = 16 generators. The generators are shown in red 
circles, the vertices of the Voronoi tessellation obtained from 
the generators in open black circles, and the edges joining the 
vertices in solid lines. (B) A centroidal Voronoi tessellation 
(CVT) obtained by applying Lloyd’s algorithm to iteratively 
update the generator positions shown in (A). The dark red 
circles show the final positions of the generators, while the 
initial position is shown by the light pink circles. The path 
traced by the generators during the iterations is shown by the 
red lines. Note that some paths cross the periodic boundaries. 
The open black circles show the vertices of the final CVT. 

To create physically realistic graphene GBs with dis¬ 
location density as close as possible to the geometrically 
required density, we propose an algorithm based on the 
centroidal Voronoi tessellation (GVT). Before describing 


the algorithm in detail, we give an intuitive explanation. 
A triangular lattice can be associated to the graphene 
lattice via a Voronoi construction (also known as the 
Dirichlet or Weigner-Seitz construction, or the dual con¬ 
struction), and vice-versa. For example, in Fig. 4B the 
graphene lattice (open black circles) forms the vertices of 
the Voronoi cells of the triangular lattice (dark red cir¬ 
cles), and vice-versa. As discussed earlier, it is difficult 
to anneal a graphene GB by using classical Monte Garlo 
methods. If however, one could find a method to anneal 
the associated triangular lattice, then the graphene GB 
could be easily recovered from the well annealed trian¬ 
gular lattice by applying the Voronoi construction. No¬ 
tice that annealing the triangular lattice by using Monte 
Garlo or MD will be almost as difficult as annealing the 
original graphene lattice with similar techniques. The 
interesting part of our algorithm uses the GVT to effi¬ 
ciently anneal the triangular lattice, and a well annealed 
graphene GB is recovered from it via a Voronoi construc¬ 
tion. 

We give a brief introduction to GVTs; details can be 
found in any number of references including Refs. [44, 45] . 
Let X = be a set of n points in a compact con¬ 

nected region C (the generalization to is anal¬ 
ogous). The points will be called the generators of the 
Voronoi tessellation. The Voronoi region Qi correspond¬ 
ing to the generator is defined as the set of all points 
that are closer (or equidistant) to it than to any other 
generator, i.e., = {x G O | ||x —x^H < ||x —Xj||, Vj 7 ^ 

i}, where || • || is the usual Euclidean norm. We denote the 
set of the vertices of the Voronoi regions by v^. Fig. 4A 
shows an example of a 2D periodic Voronoi tessellation 
with n = 16 generators. Glearly, the centroid of the re¬ 
gion is in general distinct from its generator x^. If 
we demand that the generators be arranged so that the 
centroids of the resulting regions coincide with their gen¬ 
erators, then we get a GVT. A GVT can also be described 
in terms of a variational problem [45] . It has been noted 
that the generators of a GVT are local or global mini- 
mizers of the following energy function 

n p 

HcyT(X) = V / ||x-Xj||dx. (2) 

i=i 

Fig. 4B shows an example of a 2D-periodic centroidal 
Voronoi tessellation with n = 16 generators. Note that 
the vertices of the Voronoi regions form a graphene-like 
hexagonal lattice, while the generators form a triangular 
lattice. In fact, it is a general property of GVTs that 
they tend to generate a tessellation with regular hexago¬ 
nal regions of equal size [44, 45]. The tessellation shown 
in Fig. 4B is obtained by starting from the configuration 
of generators shown in 4A and moving them according 
to Lloyd’s algorithm so as to minimize the energy func¬ 
tion Hc'yT(X) [44, 45]. The path traced by each gener¬ 
ator under the action of this algorithm is shown by the 
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red lines in Fig. 4B. We choose the aspect ratio of the 
domain such that it is possible to tile it with 16 regular 
hexagons. The tessellation is 2D-periodic because we im¬ 
plement 2D-periodic boundary conditions in our metric 
II • ||. A perfect tessellation with equal regular hexagons 
does not always exist, and neither does Lloyd’s algorithm 
converge to it from every possible initial condition even if 
it exists. For example, if we take n = 15 generators, then 
a perfect tessellation is impossible. In such cases, CVTs 
try to minimize the deviation from perfect hexagons, and 
on most occasions find a tessellation containing suitable 
pentagons-heptagons defects. 

The CVT based algorithm for generating graphene 
GBs is as follows. Given a GB the goal is to decide 
the positions of carbon atoms so that on each side of 
the GB the graphene sheet has the desired orientation, 
while the GB is comprised mostly of pentagon-heptagon 
dislocations (or undetected hexagons). To achieve this 
goal, the algorithm first generates a triangular lattice 
dual to the graphene lattice with suitable orientation on 
each side of the GB. Fig. 5 A shows an example of this 
construction. At this point, the Voronoi regions (with 
the triangular lattice points as generators) contain suit¬ 
ably aligned hexagons away from the GB, but near the 
boundary the structure is not composed of well-aligned 
pentagon-heptagon pairs, as shown in the figure. The 
generators (triangular lattice points) close to the grain 
boundaries are then relaxed by using Lloyd’s algorithm 
to obtain a GVT while keeping the points that are suf¬ 
ficiently far away from the boundaries fixed. The fixed 
points are shown by black circles, while the points that 
are relaxed by Lloyd’s algorithm are shown in red circles 
in Fig. 5. After the relaxation has converged, we obtain a 
2D-periodic Voronoi tessellation for the entire lattice, and 
obtain a graphene GB by placing a carbon atom at each 
vertex of the tessellation. Fig. 5B shows the graphene GB 
corresponding to the grain structure of Fig. 5 A obtained 
after this relaxation. We see that the grain interiors are 
defect free and have the desired orientations, while the 
GB is mediated by well-aligned pentagon-heptagon pairs. 
Finally, the atomic positions can be fine tuned by using 
the congugate gradients method and an atomistic poten¬ 
tial; we use the AIREBO potential in this study [46]. The 
graphene GB obtained after this fine tuning is shown in 
Fig. 5C. This final step only leads to small changes in the 
atomic positions, and does not entail any larger topolog¬ 
ical rearrangements of rings and defects. The numerical 
implementation is efficient, and we are able to obtain a 
well annealed GBs that are hundreds of nanometers long 
in a matter of minutes on a laptop. 

We have compared the GB structures generated with 
the proposed algorithm with the other widely used meth¬ 
ods of generating GBs. One popular technique is to paste 
together two half crystals of the required orientations and 
anneal the system by running molecular dynamics at an 
elevated temperature [26]. We use this technique, where 



FIG. 5. (A) A triangular lattice with the appropriate ori¬ 

entation is generated on each side of the GB. The triangu¬ 
lar lattice points within 10 A of the GB shown in red dots, 
while those further away are shown in black dots. Vertices 
of the Voronoi tessellation of these points make graphene-like 
hexagons in the grain interiors, but not near the boundaries, 
as seen in the figure. (B) The points within 10 A of the GB 
are relaxed using Lloyd’s algorithm in order to obtain a CVT. 
The vertices of the Voronoi tessellation of the relaxed lattice 
now comprise of hexagons and pent agon-heptagon pairs near 
the GB. The defects (pentagon-heptagon pairs) are colored for 
clarity. (C) The GB configuration obtained by placing a car¬ 
bon atom at each vertex of Voronoi regions is relaxed further 
by using the AIREBO potential and the conjugate gradient 
algorithm. The triangular lattice points are no longer shown 
for clarity. 

we heat the GB from 10 K to 3000 K, and then cool it 
back to 10 K in a span on 100 ps. The final configuration 
is then relaxed by using the conjugate gradient method. 
During this entire process a 10 A strip of atoms on the 
left and right edges of the system are constrained to their 
ideal crystalline positions. The net width of the system 
excluding the constrained atoms is 50 A. Eig. 6 shows 
the comparison of two GB structures obtained with this 
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method to those obtained by the CVT based method pro¬ 
posed here. The GBs have identical number of atoms, 
and identical atomic positions away from the boundary. 
We evaluate the GB energy per-unit length, defined as 
7 = {Etotai - natoms * -E'buik)/^^, where Etotai is the 
net potential energy of the unconstrained atoms, and 
^buik = —7.807 eV is the ground state energy per atom 
in graphene according to the AIREBO potential. It is 
evident that for the examples shown in Fig. 6 proposed 
method outperforms the method of annealing as it gen¬ 
erates GBs with lower energies. We have tested several 
hundred GBs, and the proposed method always performs 
better than the method of annealing. 


Proposed Algorithm 


Annealing 


Gm = 3.14° 

0L = 0.00° 


Gm = 38.21° 
Gl = 16.43° 



1.12 


0.80 


mization, we note that the energy landscape of the GVT 
Hamiltonian HcvtOQ is in some sense more favorable 
than that of the conventional atomistic potential based 
Hamiltonian. While just like the atomistic potentials, 
the GVT Hamiltonian can have several local minima, 
it seems that unlike the atomistic potentials, all its lo¬ 
cal minima represent low energy configurations of the 
polycrystalline graphene sheet. In a perfect tessellation, 
each generator contributes two vertices, thus removing 
or adding a generator is analogous to creating a vacancy 
pair or an adatom pair. This is a very desirable prop¬ 
erty, since it ensures that isolated vacancies or adatoms 
never appear, as these are high energy defects [47]. All 
the vacancy pair and adatom pair defects generated by 
removing and adding generators are shown in the supple¬ 
mental material and correspond to low energy configura¬ 
tions of vacancy and adatom pairs. Thus, the algorithm 
is able to produce realistic grain boundaries as well as 
point defects. 

The GVT Hamiltonian is oblivious of all the nuanced 
and complicated interactions between carbon atoms, as it 
takes a geometric view of the problem. This is a strength 
and a weakness of this approach. Its strength is clearly 
demonstrated in the high quality structures that it can 
generate at modest numerical cost. Its weakness would 
be that it is hard, if not impossible, to modify this ap¬ 
proach to include, say, the effect of chemical interactions 
with hydrogen (or another element) on structure of the 
GB. However, since the structure and properties of pure 
graphene GBs and similar two-dimensional materials are 
of such wide interest, we think that the proposed method 
has broad merit. Finally, it should be noted that the pri¬ 
mary role of the GVT algorithm is to relax the triangular 
lattice. As mentioned before, it is not feasible to simply 
use a LJ potential (or another pair potential, or hard 
spheres etc.) to relax this triangular lattice and simplify 
this algorithm. Thus, the unique properties of the GVT 
truly offer an advantage over pair potentials and Monte 
Garlo based methods in this case. 


Results and Discussion 


FIG. 6. A comparison of the GB structures and energies 
7 (in eV/A) generated by the proposed algorithm and the 
widely used technique of annealing. For both the examples 
shown here the GB generated with the proposed algorithm 
has significantly lower energy, and is thus closer to the true 
ground state than those produced with the annealing process. 
This behavior is generic, and in all the cases that we have 
evaluated the proposed method almost always results in lower 
GB energy. 

To understand why this boundary generation algo¬ 
rithm outperforms the traditional method of annealing 
or grand canonical Monte Garlo or simple energy mini- 


The full library of our measured experimental GB 
structures is available at the experimental structures 
archive. The full library of our simulated GB structures 
is available at the computed structures archive. 

Fig. 7 shows 17 examples of experimentally measured 
GB structures, ranging from low to high boundary dis¬ 
orientations II^mII- These misorientations are calculated 
from the best-fit lattices of the two grains, with an es¬ 
timated error of approximately 0.5°. The low angle 
boundaries are composed of isolated pentagon-heptagon 
pairs, while the higher angle boundaries are composed of 
connected pentagon-heptagon pairs. Each experimental 
boundary is paired with a matching example taken from 
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FIG. 7. Examples of experimentally measured graphene GBs (above) compared to similar numerically simulated boundary 
structures (below), sorted by disorientation angle || (^m +30° mod 60°) — 30°||. Boundary structures range from isolated 
dislocations, periodic arrays of separated dislocations, continuous high-symmetry boundaries to serpentine boundaries with 
large amounts of excess dislocation content. 


the generated boundary library, with either an identical 
or very similar structure. The close agreement demon¬ 
strates the efficacy of our boundary generation algorithm. 


Simulated Grain Boundary Structures 

A small subset of the numerically simulated GB struc¬ 
tures are plotted in Fig. 8 A. The 5-member pentagon 
rings are colored in red, while the 7-member heptagon 
rings are colored in blue. Each pentagon-heptagon 
pair sharing a C-C bond represents a (1,0) dislocation 
core with the smallest possible Burgers vector, while a 
pentagon-heptagon pair connected by a C-C represents 
a (1,1) dislocation core with the next-smallest Burgers 
vector [11]. Separating the pentagon and hepagon by a 
single 6-member hexagon ring leads to a (2,0) dislocation 
with a Burgers vector with twice the magnitude of the 
(1,0) dislocations. Fig. 8B shows the atomic structure of 
these dislocations graphically. 


Fig. 8B also shows another commonly observed dis¬ 
location structure; pairs of (1,0) and (0,1) dislocations. 
These dislocation pairs have the same Burgers circuit as 
the (1,1) dislocation. These dislocation pairs are typi¬ 
cally much lower energy than (1,1) dislocations and are 
commonly observed in the range 21.8° < < 60°, 

and are especially prevalent in for misorientation angles 
Om > 32.2° [11]. 


Structural Properties of Simulated Boundaries 

We have used automated analysis routines to extract 
the dislocation content and structural properties from all 
of the lowest-energy simulated boundaries for each value 
of Om and 6 >l. Figs. 9A-F plot the dislocation content of 
the simulated boundaries. The most common boundary 
structures by a large margin (over 98%) are the (1,0) 
and paired (1,0)+(0,1) dislocations. Figs. 9B and C show 
that the pairing arrangement is much more common for 










0Line = O° 10° 20° 30° 40° 50° 60° 



FIG. 8. (A) Examples of the GBs calculated with our bound¬ 
ary generation algorithm for various misorientation angles 
and boundary line angles ^l- (B) Dislocation structures 
present in low-energy graphene GBs. 


boundaries with Om > 30°, although many pairs are also 
present for the most asymmetric boundaries (high 6 >l). 

Figs. 9D and E plot the density of (1,1) and (2,0) dis¬ 
locations, both of which are almost entirely present only 
in boundaries with high disorientations, 25° < Om < 35°. 
The peak density of (1,1) and (2,0) dislocations is approx¬ 
imately 10 and 60 times lower than the (1,0) dislocation 
density respectively, shown in Fig. 9F. (1,1) dislocations 
are slightly more prevalent at higher Oj, values, while (2,0) 


dislocations have higher density at lower 6 >l values. 

We have also analyzed the three-dimensional orienta¬ 
tion densities of the dislocation dipole vectors, defined as 
the vector from the center of each heptagon to its asso¬ 
ciated pentagon. The 2D probability distributions of all 
dislocations (equally weighted for each calculated bound¬ 
ary) over in-plane and out-of-plane dipole tilt vectors are 
plotted in Figs. 9G, H and I for (1,0), (1,1) and (2,0) 
dislocations respectively. The (1,0) dislocations tend to 
align along the boundary line, with two large clusters 
visible in Fig. 9G; the left cluster is formed from the 
lower boundaries while the cluster to the right con¬ 
tains more high Om boundaries. These right-side (1,0)- 
type dislocations tend to be slightly tilted away from 
the boundary line vector and are often paired with a 
(0,l)-type dislocation, giving a longer tail towards lower 
values in this cluster. A third, very dim cluster is 
visible at approximately 60° in-plane tilt values. All 
three of these clusters are centered on 0° out-of-plane 
tilt, with the distributions decreasing quickly at higher 
and lower out-of-plane tilt values. All three clusters have 
a range of approximately ±30° for the out-of-plane tilt. 
By contrast. Fig. 9H shows that the (1,1)-type disloca¬ 
tions have a strongly bimodal probability distribution for 
both in-plane and out-of-plane tilts, occur primarily at 
in-plane tilts of ?^20° and ^170° and out-of-plane tilts of 
±13.5°. The (2,0)-type dislocation orientations are plot¬ 
ted in Fig. 91, showing maxima at an in-plane tilt of 0° 
and out-of-plane tilts of ±12.5°. 

The roughness of all simulated boundaries was esti¬ 
mated by connecting all boundary pentagons and hep¬ 
tagons sequentially, and measuring the root-mean-square 
(RMS) displacement of this distorted boundary line, both 
in-plane and out-of-plane of the graphene sheet. The in¬ 
plane RMS roughness is plotted in Fig. 9J as a function of 
the boundary angles. The mean and standard deviation 
of the roughness averaged over all 6 >l values is plotted in 
Fig. 9K. The in-plane boundary roughness is largest at 
low values, decreasing from approximately 2.5A to 2A 
with increasing ^m- The in-plane roughness of the sym¬ 
metric tilt boundaries are also plotted in Fig. 9. Most 
symmetric tilt boundaries have lower roughness than the 
average of all boundaries, approximately lA, except for 
a small number of boundaries spiking at RMS rough¬ 
ness values of 3A. The out-of-plane RMS roughness of all 
boundaries is plotted in Figs. 9L and M. These roughness 
values are much more uniform that the in-plane rough¬ 
ness; between Om = 0° and 15° the roughness decreases 
from 2 . 4 A to 1.3A. From = 15° and 45° the roughness 
is almost constant at 1.3A. Finally, between = 45° 
and 60° the roughness increases from 1.3A to 2.4A. The 
symmetric boundaries do not show any deviation in out- 
of-plane roughness compared to the average values. The 
boundaries with larger disorientations have higher dislo¬ 
cation densities; this allows their overlapping strain fields 
to more easily cancel out and therefore lead to lower out- 
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FIG. 9. Dislocation density probability distributions for (A) (1,0), (B) paired (1,0)-(0,1), (D) (1,1) and (E) (2,0) dislocations. 
(F) Average dislocation density as a function of the misorientation angle for (l,0)-type dislocations showing fractions of 
primary (1,0) and secondary (0,1) orientations. (F) Logarithmic probability of average dislocation densities for all types. 
Dislocation dipole vector orientation probability distributions for (G) (1,0), (H) (1,1) and (I) (2,0) dislocations. Boundary 
line RMS roughness for (J)-(K) in-plane and (L)-(M) out-of-plane displacements. Average roughness is plotted as red line, 
while the standard deviation is plotted as a pink boundary. 


of-plane roughness. 


Structural Properties of Experimental Boundaries 

The physical properties and structure of the exper¬ 
imental boundaries were also characterized with auto¬ 
mated routines. These results are shown in Fig. 10. 
Some examples of symmetric tilt boundaries with struc¬ 
tures similar to those plotted in Fig. 7 are depicted in 
Fig. 10 A. The experimental results plotted in the rest 
of Fig. 10 are 160 boundaries estimated to be unique 
structures. We took this step to try to minimize dou¬ 
ble counting of boundary structure datapoints. We ob¬ 
serve that most of the boundary structures we measured 
fall at high disorientation angles, i.e. boundaries close 
to = 30°. This phenomenon is due to topological 
effects; low angle boundaries have rougher surfaces and 
long range out-of-plane distortions [12]. This topology in 
turn attracts carbon contamination due to surface charg¬ 
ing, which obscures the boundary. The imaging process 
is therefore biased towards flatter boundaries, i.e. those 
closer to Om = 30°. 

Figs. lOB, C and D shows the measured densities as a 
function of misorientation angle for (1,0), (1,1), and 
(2,0) dislocations respectively. These figures also show 
the (1,0) dislocation densities of the 6 symmetric bound¬ 


aries plotted in Fig. 10, and the dislocation densities of 
all three types predicted from the simulated boundary re¬ 
laxations in Fig. 9F. The experiments are in good agree¬ 
ment with both of these sets of predictions. All of the 6 
symmetric boundaries shown in Fig. 10 A have a nearby 
experimental example. However at misorientation angles 
in the range 25° < < 35°, in the highest density 

region of the experimental boundaries, the average pre¬ 
dictions of the relaxed and constructed boundaries are 
much closer to the majority of experimental dislocation 
densities. The simulated boundaries also predict a small 
concentration of (1,1)- and (2,0)-type dislocations in the 
range 25° < Om < 35°, both of which are observed in 
the experimental measurements shown in Figs. IOC and 
D. The average dislocation concentrations of the experi¬ 
ments are very close to the simulations, with the excep¬ 
tion of a single (1,1) dislocation observed in an experi¬ 
ment at Om = 14°. 

Because the experimental boundaries are measured as 
a 2D projection, the out-of-plane distortions cannot be 
directly measured. However, these distortions are typi¬ 
cally accompanied by large deviations in the local pro¬ 
jected atomic positions. We have therefore measured 
the average strains and local rotations of all boundaries, 
plotting in Figs. lOE and F, normalized by the dislo¬ 
cation density. All of the average strain metrics have 
approximately the same trend; they decrease as the mis- 
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A 137 119 17 113 131 191 



FIG. 10. (A) Examples of special symmetric bound¬ 

aries with the smallest repeat lengths. Experimental dislo¬ 
cation density for (B) (1,0), (C) (1,1) and (D) (2,0) dislo¬ 
cations respectively. (E) RMS atomic strain of experimental 
boundaries, parallel and perpendicular to boundary line. (F) 
RMS atomic rotations of experimental boundaries. Trendlines 
shown for experimental strain and rotation. (G) In-plane 
RMS roughness of experimental boundaries and the bound¬ 
aries in (A). 


orientation increases, and then the strain increases past 
Om = 30°. This is qualitatively in agreement with the 
predictions of out-of-plane roughness trends shown in 
Fig. 9M. The RMS strain perpendicular to the grain 
boundaries is large than the parallel strain for virtu¬ 
ally all boundaries, shown in Fig. lOE. This is because 
dislocation dipoles are typically aligned along the grain 
boundaries, which allows the adjacent dislocation strain 
fields to partially cancel out. 

The boundary RMS roughness for all experimental 
boundaries is plotted in Fig. lOF. Five of the six sym¬ 
metric boundaries plotted in Fig. 10 A predict boundary 
roughness values that are very close matches to the exper¬ 
iments. The largest concentration of boundaries near the 
center of the plot reach a minimum roughness at a mis- 
orientation value closer to the E13 value of = 32.2°. 
The entire cluster of values has an average in-plane RMS 
roughness of approximately 2.5 A, in good agreement with 
the synthetic boundary prediction of 2Agiven in Fig. 9K. 

Matching Experimental and Simulated Structures 

The vast majority of experimentally measured GBs 
and all of the numerically simulated GBs in this study 
can be constructed by mixing the dislocations structures 
shown in Fig. 8B. Symmetric boundaries with < 21.8° 
contain aligned (1,0) dislocations, while non-symmetric 
boundaries and boundaries with Om > 21.8° are typically 
composed of a mixture of (1,0) and (0,1) dislocations [11]. 
The ratio of the number of the two most common orienta¬ 
tions for (1,0) dislocations could be used to estimate the 
boundary line angle but for many of the experimental 
images the boundary length is too short (not enough ob¬ 
served dislocations) for an accurate measurement of the 
line angle. 

As predicted, the low angle boundaries consist of iso¬ 
lated (1,0) dislocations. At low misorientation angles, all 
of the (1,0) dislocations have the same orientation, while 
at high misorientation angles (e.g. Om = 48°) the struc¬ 
ture consists of (1,0) and (0,1) dislocation pairs. Three 
of the plotted examples in Fig. 7, Om = 20.8°, 21.3, 
and 22.2, have structures very similar to the E7 special 
boundary [11, 48, 49]. The Om = 32.6° boundary is a 
nearly perfect example of the 1113 special boundary [11]. 

The high angle boundaries in Fig. 7 are formed from 
continuous or near-continuous dislocation groups with al¬ 
ternating 5- and 7-member rings. This structure is ex¬ 
pected, since a more equal local density of pentagons 
and heptagons leads to a lower net disclination content 
and thus less 3D topological variation, and more sta¬ 
ble structures [12]. The longest high angle boundaries 
in Fig. 7 {Om = 30.4°,30.2°) show an interesting devi¬ 
ation from the flat boundary structures; they form ser¬ 
pentine structures similar to some literature predictions 
and observations [20, 31, 50] as well as our previous ex- 
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periments [3, 35]. Both of these boundaries exhibit two 
relatively sharp 30° deviations from a flat boundary line, 
where the arm-chair and zig-zag graphene edges of the 
two sides exchange identities. These structures likely 
originate from capillary fluctuations during the initial 
growth, but form (relatively) well-defined faceted edges 
rather than a rougher boundary. 

Overall, the boundary generation algorithm therefore 
accurately predicts structural properties of the exper¬ 
imental boundaries. The primary disagreement is the 
slightly higher boundary roughness and dislocation den¬ 
sity of the high angle boundaries near ^m- The source 
of this minor disagreement lies in the faceting exhibited 
by the experimental boundaries, such the bottom two ex¬ 
perimental structure plots in Fig. 7. Since the simulated 
boundaries were constrained to follow a single line an¬ 
gle (flat boundaries), they could not fully capture this 
effect. However, our CVT algorithm could easily be used 
to simulate such boundaries. This example shows how 
large-scale experiment and simulations couple together, 
where the different areas of agreement or disagreement 
can be used to improve structure models. 


Conclusion 


In summary, we have used semi-automated processing 
routines to characterize the structure of a very large num¬ 
ber of experimentally measured single-layer graphene 
grain boundaries and described their local atomic struc¬ 
ture as a function of misorientation angle. We have 
also introduced a new algorithm for generating realis¬ 
tic graphene grain boundaries that produces structures 
in excellent agreement with the experimental boundary 
structures. We have used a combination of our algorithm 
and molecular dynamics to generate and relax graphene 
grain boundary structures covering the entire orientation 
parameter space for single-later graphene boundaries. 
The structure and physical properties of these simulated 
boundaries were analyzed as a function of misorienta¬ 
tion and boundary line angles. A detailed comparison of 
the experimental and simulated boundaries demonstrates 
that our structural models have high predictive power for 
the experimental structures. In a forthcoming study, we 
will analyze the energetics of both experimental and sim¬ 
ulated grain boundaries. Finally, all experimental and 
simulated structures are made available on the internet. 
We hope that this paradigm of computer-assisted anal¬ 
ysis of a statistically relevant number of structures and 
availability of all measured data becomes standard for 
the study of atomic-resolution structures. 
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